clear; close all; clc;

%% 参数设置
% memcapacitor 参数
alpha = 2.1;              % 控制参数 α
beta  = 2;                % 控制参数 β

% 激励信号 q 的幅值
q_amp = 2;

% 选择不同的频率
F_list = [0.1, 1, 20];

% 仿真步长
dt = 0.001;

%% 绘图
figure;

% 子图1: memcapacitor 的频率相关钳制滞回特性（Fig1(a)）
subplot(1,2,1);
hold on;
colors = ['b','r','g'];  % 为不同频率选择不同颜色

for idx = 1:length(F_list)
    F = F_list(idx);
    T_period = 1/F;        % 一个周期
    T_total  = 30 * T_period; % 总仿真时间（多周期确保稳态）
    
    t = 0:dt:T_total;      % 时间序列
    q = q_amp * sin(2*pi*F*t);  % 输入信号 q(t) = 2*sin(2π·F·t)
    
    % 用欧拉法求解微分方程 dσ/dt = q - σ，初始条件 sigma(1) = 0
    sigma = zeros(size(t));
    for k = 1:length(t)-1
        sigma(k+1) = sigma(k) + dt * (q(k) - sigma(k));
    end
    
    % memcapacitor 电压计算: v_C = (alpha + beta*sigma) * q
    vC = (alpha + beta * sigma) .* q;
    
    % 取最后一周期的数据进行绘图，以避免初始瞬态
    ind_start = find(t >= (T_total - T_period), 1, 'first');
    plot(q(ind_start:end), vC(ind_start:end), 'Color', colors(idx), 'LineWidth', 1.5);
end

xlabel('q');
ylabel('v_C');
title('memcapacitor 频率滞回曲线');
legend('F = 0.1','F = 1','F = 20','Location','Best');
grid on;
hold off;

% 子图2: N型LAM memristor 的直流 V-I 曲线（Fig1(b)）
subplot(1,2,2);
gamma = 2;   % memristor 参数
% ϕ 取值范围，选取足够区间保证包含局部负斜率部分
phi = linspace(-1, 1, 1000);

% 根据模型： V = (ϕ - tanh(ϕ))/gamma, I = (ϕ^2 - 0.5) * V
V_dc = (phi - tanh(phi)) / gamma;
I_dc = (phi.^2 - 0.5) .* V_dc;

plot(V_dc, I_dc, 'm', 'LineWidth', 1.5);
xlabel('V');
ylabel('I');
title('N型LAM memristor 直流 V-I 曲线');
grid on;
